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Abstract 

We show that correlated dynamics and long time memory persist in self-organized 
criticality (SOC) systems even when forced away from the defined critical point 
that exists at vanishing drive strength. These temporal correlations are found for 
all levels of external forcing as long as the system is not overdriven. They arise 
from the same physical mechanism that produces the temporal correlations found 
at the vanishing drive limit, namely the memory of past events stored in the system 
profile. The existence of these correlations contradicts the notion that a SOC time 
series is simply a random superposition of events with sizes distributed as a power 
law, as has been suggested by previous studies. 
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1 Introduction 



Are events in a self-organized criticality (SOC) [1,2] system random or cor- 
related? This is a still ongoing and important topic of discussion in the SOC 
community that has been argued from both sides to propose or invalidate 
practical tests for SOC behavior [3,4,5,6,7,8,9,10,11]. One specific point of 
this dialogue concerns the distribution of the amount of time between the 
initiation of successive events (waiting, laminar or first return times) in a sys- 
tem suspected of exhibiting SOC behaviour. At issue is whether or not the 
probability density function (PDF) of this measure should always scale expo- 
nentially or as a power law for a SOC system. Three recent examples of this 
arc drawn from space plasma physics [3], laboratory confined fusion plasma 
physics [4] and earthquake geophysics [9]. The authors of all of these studies 
first assume that events in a SOC system must be uncorrelated. They then 
report observations of waiting times distributed as power laws and draw sep- 
arate but similar conclusions that the systems' (space plasma, fusion plasma, 
earthquakes) dynamics cannot be consistent with SOC. 

The root of this controversy is to be found in the misconceptions that "strong 
correlations between successive [events is] at variance with the SOC model" [3] 
and that an event "can occur randomly anywhere at any time and. . . cannot 
'know' how large it will become" [7,9]. Based on these statements, the idea 
of exponentially distributed waiting times could be easily derived [3]. But it 
has been known for some time that these notions simply do not hold in a 
SOC system [12,5]. They were probably encouraged, in the early days of SOC, 
by analytical calculations [13,14] that appeared to reproduce quite well the 
power spectrum of the avalanche activity observed in the Bak-Tang-Wiesenfeld 
(BTW) sandpile [1] (and later in the directed running sandpile [15]) when 
modeling it as a random superposition of events with power law distributed 
sizes. Since then, a state of confusion that still persists appears to surround 
this issue, at least among those looking for SOC features in physical systems, 
with quite a few papers accepting [3,4,9,11] and others denying [5,8,10,16] the 
validity of these notions. 

In this paper we show explicitly that features exist in the power spectra of 
SOC systems that signify dynamical correlations and that, therefore, any as- 
sumption or test based on the absence of temporal correlations is meaningless. 
By 'dynamical' we mean that temporal correlations are caused by the system 
dynamics, not by a hypothetical non-randomness in the external drive [5]. To 
study them, we have examined the avalanche activity time series of the same 
one dimensional running sandpile model analyzed in Ref. [15] as our prototyp- 
ical SOC model. The benefit of using this model is that it allows us to analyze 
systems with a finite drive strength instead of just the vanishing drive limit in 
which the absorbing phase transition responsible for the appearance of crit- 
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icality and universal behaviour was identified [17,18,19,20]. Indeed, in many 
practical applications a finite drive is more appropriate. This is the case, for 
instance, in most magnetized plasmas, either inside the Sun [21] or confined 
in an earth-based tokamak [22]. 

Taking advantage of increased computer power over that available at the time 
of Ref. [15] we explore timescales up to four decades longer than those used 
in that study for varying driving rates and system sizes. Because of this, dy- 
namical correlations not previously detected in the power spectrum of the 
avalanche activity time series (associated with f"^ power law regions of ex- 
ponent < (3 < 1 ) are now apparent. The exponent is found to increase with 
the drive strength, reaching (3 = 1 only at the higher drives. We show that 
this is not a signature of the increasing dominance of avalanche overlapping 
(as stated previously in Ref. [15]), but just a consequence of the fact that 
the timescales where the correlations are present are also displaced towards 
higher frequencies at stronger drives. Therefore, avalanche overlapping is not 
the cause of the 1// region seen in the directed running sandpile. 

The changing value of /? might seem to suggest that the character of the 
correlations is changing with the drive as we move away from the critical point. 
However, when analyzing the same activity signals with rescaled range {R/S) 
analysis [23,24] to estimate the Hurst exponent H, a very different conclusion is 
drawn. Dynamical correlations are associated with H > 0.5. We find a constant 
Hurst exponent H ~ 0.8 for all values of driving rate over the same timescales 
described above, where /3 > in the power spectrum. H remains constant 
as f3 changes from /3 ~ 0.4 at the lowest driving rate studied here to /5 ~ 1 
at the highest driving rates. Interestingly, in the vanishing drive limit [25] we 
find the same value of H. In this natural limit of the running sandpile model, 
avalanches do not overlap, by definition, because the external drive is turned 
off as soon as an event is initiated. The origin of the correlations must thus 
be found, then, in a mechanism other than avalanche overlap. Specifically, the 
origin of temporal correlations in the sandpile model is the shape of the system 
profiles, carved by past avalanches, where the memory of the system history 
is stored. 

In this report, we do not set out to prove that temporal correlations exist in a 
SOC system. In the words of a helpful referee, "there is no question whether 
correlations exist in self-organized or any other critical systems. They do, ba- 
sically by definition." Such correlations can be proved through the existence of 
diverging correlation lengths and times [26] . Given that temporal correlations 
exist in SOC systems, our goals are to lay to rest the controversy mentioned 
in the first paragraph and to address the questions: How does one identify 
temporal correlations in a real system? What measures should one use? What 
do those measures say? What mechanism(s) produce(s) the correlations? We 
use the power spectrum and R/S analysis in this study but do not claim that 
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they are the only possible applicable measures. Multiple measures should be 
used for completely thorough studies of dynamical systems. 

The paper is organized as follows: in Sec. 2 the basics of the running sand- 
pile and analysis techniques are reviewed. Numerical results obtained for the 
avalanche activity time series arc described in Sec. 3 and discussed in Sec. 4. 
Finally, some conclusions are drawn in Sec. 6. 



2 Model and measures 

2.1 Model 

The running sandpile can be described by four parameters: L, Pq, Zcrit and 
A^f. Consider a single column of L cells. Each cell contains an integer number 
of "sand grains" ; this number is the height of the cell. Sand is added to each 
cell by a random "rain" from above. That is, at each time step for each cell, 
there is a probabihty < Pq < 1 that a grain of sand will be added to it. 
The average input current into the entire system is Jin = PqL grains per time 
step. The local gradient is the difference in height between two neighboring 
cells. If this local gradient exceeds a critical gradient Zcrit then an avalanche 
occurs. An avalanche stabilizes the local gradient by transferring A^f grains of 
sand from the higher cell to the lower in a single event called a flip [25] or 
a toppling (the more standard term). This avalanche can make .^n+i and/or 
Zn_i unstable at the next time step so that the avalanche spreads to other 
cells. In this way, spatially and temporally extended events in a system can 
occur. This sandpile can be extended to two or more dimensions [2] but we 
will only discuss the one dimensional case. For this study, we want to explore 
the effects of changing system size (L) and external driving rate (probability 
Pq) so we have arbitrarily fixed the parameters Zcrit = 8 and A^f = 3. In 
Ref. [22] it was shown that the system dynamics does not qualitatively change 
for different values of Z^rit and iVf, as long as 1 < iVf < Zcrit/2. 

Note that time is well-defined in this running sandpile model as opposed to the 
vanishing drive model (the original BTW model, for instance) where external 
forcing (sand addition) is suspended as long as there is activity in the system. 
In this running sandpile, sand can be added at each cell at every time step 
with a probability Pq even if an avalanche is active. When we increase the 
level of external forcing by increasing Pq, we increase the probabilty of a 
grain of sand being added to each cell and therefore increase the probabilty 
of an avalanche. But we do not change the time step — the clock continues 
to tick regularly regardless of system activity. Computationally, this requires 
simultaneous random adding of grains (forcing) and checking for instability 
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(relaxation) within the same programming loop. A natural time scale for the 
system, then, can be thought of as Pq^- 

The time series that we analyze is the instantaneous total avalanche activity 
or, in short, the activity pa, which is defined as the total number of unstable 
sites (where Z > Z^rw) at each time step. An unstable cell transports A^f 
grains of sand to the next cell. The activity at each time step can be thought 
of as the instantaneous (potential) energy dissipation in the system or the 
total instantaneous flux through it. Avalanches in a Pa time series appear as 
structures — sequential nonzero activity separated by periods of inactivity, or 
quiet times. We analyze activity data for times after steady state has been 
reached, which is when the long time average number of topplings approaches 
a constant value. 

Two temporal measures from the activity series are discussed in this paper: 
durations and quiet times. A duration is the number of time steps during 
which the system is active; we refer to this as an avalanche. Note that for the 
running sandpile that we use in this study, an avalanche defined in this way 
could actually comprise two or more distinct avalanches that happen to be 
active at the same time in the system. At high drive, this tends to be the case 
but at low drive individual avalanches can be distinguished through analysis 
of the pulse shape in the activity time series. A quiet time is defined as the 
length of an uninterrupted sequence of zero topplings (inactivity) between 
two avalanches. It is the length of time that the system is in an absorbing 
state. In the randomly driven sandpile, probability density functions (PDFs) 
of unadulterated avalanche durations scale as a power law in the self-similar 
range [27] while the PDFs of quiet times scale exponentially [6,5]. 

An important part of this study is to quantify the effect of changing the drive 
strength and system size on correlations. First, it is important to note that 
for critical dynamics to exist, the system must not only be underdriven (i.e., 
the flux in, Jjn = -Pq-^) must never exceed the maximum possible flux out of 
the bottom cell, Jqut = Af)? b^it it must also satisfy Jin < so that 

a cell can alternate between stable and unstable until the avalanche ends or 
washes past the cell [22]. Second, note that the difference between driving 
close/far from the condition Jin < Af/2 is to increase/decrease the degree 
of avalanche overlapping. That is, if Jin ^ A^f/2 then there will be very few 
cases of simultaneous avalanches and hence very little avalanche overlapping. 
In practise, Jin can be set as small as desired to effectively eliminate all over- 
lapping, thus allowing identification and analysis of single avalanches. In this 
limit, avalanche dynamics and statistics are the same as the vanishing drive 
(non- running) sandpile of [25]. 

We want to quantify the amount of temporal event overlapping and we want to 
compare systems of different sizes (L). If the average avalanche duration d (~ 
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Fig. 1. Partial activity time series from running sandpile model at steady state for 
(a) low and (b) high drive. At low drive {PqL'^ = 0.2), there are approximately 30 
distinct and separate events in ~ 10^ time steps while at high drive {PqL? = 100) 
there are far more events and their identification is made difficult by almost constant 
overlapping. Note the trapezoid shape of a single, distinct avalanche in the magnified 
inset of the low drive case. A superpulse is defined for high drive as the sequence of 
activity between consecutive quiet times. 

L) is greater than the average trigger time Tt (~ Ni/P^L) [5] then avalanches 
will overlap and therefore the condition for overlapping in time is d > Tt or 
PqL^ > N{. (This is in contrast to the limit for spatial overlap previously given 
in Eq. (4) in [15].) For a fixed Nf, this condition encompasses both size and 
drive strength so we define the effective driving rate as Je = PqL^- At steady 
state, of course, flux into and out of the system are still equal, Jin = PqL. 
Je is simply a parameter defined to compare systems of different sizes that 
have proportional levels of overlapping. In what follows, low drive will refer 
to systems where overlapping of events is negligible (i.e., Je -C iVf), and high 
drive will refer to cases when overlapping becomes significant (i.e., Je ^ Nf) 
(see Fig. 1). This will be discussed further in Section 5. 



2.2 Measures 



Though a variety of measures can be used, our tools to detect correlations are 
power spectra and rescaled range {R/ S) analysis. The discrete power spectrum 
of a data set X{t) is defined as the square of the discrete Fourier transform, 
= where 

N-l 

F{f) = AT-i ^ X(i)e-^2'^(^/^)*. (1) 
t=o 

What does the power spectrum say about correlations in a series? More ap- 
propriate here, what does a power spectrum that scales as a power law f^^ 
have to say about persistence in a series? For (3 > 0, the low frequencies are 
more important in the dynamics so the series tends to look relatively smooth 
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and trends persist. Such series are said show persistence. For /3 < 0, the high 
frequency components are more important and the series is very rough. Trends 
do not persist for long and the series is called antipersistent. Truly random 
noise (white noise) has a spectrum that scales with P — Q [28] . 

The Hurst exponent, H G [0, 1] [23], can be used to measure persistence in a 
time series. A value of H > 0.5 implies persistence, H < 0.5 antipersistence 
and H = 0.5 uncorrelated. We measure H by using R/ S analysis [23,29,24]. 
The rescaled range is defined as R'{t) = R{t)/S{t), where S{t) is the stan- 
dard deviation and 

R{t) = maxPF(A;,r) — min W{k,T) (range), 

l<fe<T l<fe<r 

k 

W{k,r) = ^(^t — {X)t-) (cumulative deviation) and 
t=i 

{X),^-j^Xt (mean). 

If the rescaled range of the time series exhibits regions that scale as R'{t) ~ 
T^, then H is the Hurst exponent. 

Need to work on this paragraph. In some cases, H as estimated by R/ S 
analysis can be related to (5. For instance, for a time series of stationary 
independent identically distributed random variables (drawn from a Gaussian 
distribution), (3 = 2H — 1 [30]. If, however, the random variables are instead 
distributed with a fat tail that scales as ~ then the two parameters are 
more generally related as /3 = 2{H + 1 — A*"^) [24]. In both examples, a crucial 
underlying assumption is that the random variables that comprise the time 
series are stable and stationary. This should be kept in mind for Section 4.3 
below, where we show and discuss the fact that /3 and H are not simply related 
for the activity of the running sandpile. 



3 Results 

In this section, we describe the primary qualitative features of the power spec- 
tra and R/S analysis of activity time series obtained for different driving rates 
and system sizes. The main theme is that distinct power law regions always 
exist in both measures; we defer the interpretations of the regions to the fol- 
lowing section. 

The power spectra of the activity look very different depending upon whether 
the driving rate is low or high. Fig. 2 shows the spectra for over five orders of 
magnitude of effective driving rate, Je = PqL^, which increases from top to 



7 



: I iiiin 1 1 1 1 mill 1 I Mini 1 1 1 iiiiiii 1 i iiiin 




Fig. 2. (Color online.) Power spectra of activity time series of L = 200 sandpile 
for five orders of magnitude of effective driving rate in PqL^ € [0.002, 296]. Spectra 
have been shifted along the y axis for easier viewing. 

bottom in the figure. [The sandpile size is L = 200, but very similar results 
have been obtained for sizes between L = 50 and L = 2000.] The lowest drive 
used is PqL'^ = 0.002 and the highest is PoL"^ = 296, chosen to stay below 
the normal overdrive limit of PqL < N{/2, discussed above. Multiple distinct 
power law regions are seen in all of the spectra. The highest frequency regions 
share a common slope (~ 3.5). At low drive a very prominent hump at low 
frequency moves to higher frequency as driving rate increases. 

For clarity, selected spectra for a low, a medium and a high drive case are 
shown in Fig. 3(a) for L = 200. For the lowest drive case, six power-law 
regions can be distinguished, labeled A to F as frequency decreases. They 
are separated by breakpoints denoted by T^^, k being the appropriate region 
label. Similarly, four regions are found in the highest drive case, labeled A, 
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Fig. 3. (Color online.) (a) Power spectra and (b) R/S analysis of activity for fixed 
system size and low, medium and high driving rates. The y values of both measures 
have been shifted for easier viewing. In the spectra {R/S), six (five) regions of low 
drive and four regions of high drive are shown by solid lines. Lines are power laws; 
numbers show values of f3 in spectra and H in R/S. Lowest frequency region of 
spectra and H = 0.5 of R/S for the low drive case is not seen because of the finite 
size of the time series. Its existence is assumed based on the regions seen in the 
spectra of higher drive cases. 




Fig. 4. Cartoons of distinct regions and their breakpoints and causes of power spectra 
and R/S analysis of sandpile activity, (a) Power spectrum of low drive, (b) R/S 
analysis of low drive, (c) power spectrum of high drive and (d) R/S analysis of high 
drive, (c) is taken from Fig. 6 of [15] and the others are drawn in that spirit. 
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Fig. 5. (Color online.) (a) Inverse breakpoints of the power spectra versus break- 
points of R/ S analysis of activity for different size sandpiles and driving rate and 
(b) breakpoints of power spectra versus driving rate for L = 200 sandpile. Numbers 
shown are exponents of power law fits to the data. 

Z), E and F for reasons to be discussed in the next section. The regions and 
the labelhng convention are illustrated in the cartoons in Figs. 4(a) and (c). 
[Fig. 4(c) is based on Fig. 6 of Ref. [15].] 

Next, we will describe the R/ S analysis results. They are shown, for the same 
selected cases, in Fig. 3(b). Again, different power law regions can be distin- 
guished: five distinct regions for the lowest drive case and four for the high- 
est. They are labeled according to the same convention as the power spectra 
(Figs. 4(b) and (d)), except for the shortest timescale region, labeled instead 
as "A/B". Again, relevant breakpoints are identified and labeled. 

The labelling chosen for the spectrum and R/ S analysis breakpoints is not 
arbitrary. Breakpoints are compared in Fig. 5(a), where it is found that those 
of the two measures, found independently, agree very closely with each other. 
[The R/ S breakpoints appear at slightly longer timescales than those of the 
power spectrum, but this effect is known from comparisons with Hurst ex- 
ponents determined via the structure functions method [32,31].] We conclude 
that both measures distinguish the same dynamical regions through the iden- 
tification of different power law regions. 

For the discussion in the next section it is interesting to know how these 
breakpoints scale with driving rate (see Fig. 5(b)). Two different behaviors 
can be distinguished: the breakpoints associated with the A and B (and A/B) 
regions are independent of the drive, while all of the others scale (almost) 
linearly with it, moving to higher frequencies (or shorter time lags for the 
R/ S analysis) as the drive increases. 

To conclude the enumeration of results, we will briefly describe the effect 
that changing the sandpile size has on the signatures and then more fully 
discuss this effect in Section 5. If L is increased while PqL"^ is held constant, 
then all regions move to smaller frequencies (or longer time lags). Fig. 6. The 
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Fig. 6. Power spectra and R/S analysis of activity for low drive sandpile for 
PqI"^ = 0.2 and L = 50, 100, 200, 1000. Plots have been shifted along the y axes 
for better viewing. The solid lines are power laws plotted with the labeled expo- 
nent; they are not fits to the data. These plots show how the distinct regions move 
to longer time scales with larger system size but with the same effective driving 
rate, PqLP' . 
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Fig. 7. Power spectra and R/S analysis of low drive sandpile: original series, with 
quiet times removed and with pulses shuffled. Plots have been shifted along the 
y axes for better viewing. The solid lines are power laws plotted with the labeled 
exponent; they are not fits to the data. The spectrum and R/S of the activity time 
series for the zero drive sandpile are virtually identical to the spectrum and R/S 
of the low drive series with quiet times removed; to preserve clarity, they are not 
shown. 

values of the exponents associated with each distinct region remain essentially 
unchanged, especially in the R/ S analysis. 



4 Discussion 



The first point that can be made from Figs. 2, 3 and 6 is that long time cor- 
relations do exist in the running sandpile regardless of driving rate. Changing 
the rate at which sand is dropped only changes the timescales over which the 
correlations are seen. This is evident in the spectra of Fig. 2, where the promi- 
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nent hump, whose right side contains the (0 < /3 < 1) correlated region, 
moves from low frequency at low drive to high frequency at high drive. The 
fact that this hump exists at all is proof of correlations in the system, since 
a random superposition of pulses has a flat /° spectrum at low frequencies 
[13,14]. This is further confirmed by taking a activity time series from a low 
drive case, where events are distinct and separated, and randomly shuffling 
the events to make a new time series. The spectrum and R/S of the new 
time series show a /° and H = 0.5 region, respectively, indicating no corre- 
lations, whereas before both measures showed correlated dynamics (Fig. 7). 
Note that this simple experiment also tells us that it is the avalanche order in 
which correlations are present, since the order of the quiet times among them 
merely reflects the random character of the drive [6]. Therefore, the previously 
mentioned notions that "strong correlations between successive [events is] at 
variance with the SOC model" [3] and that an event "can occur randomly 
anywhere at any time and. . . cannot 'know' how large it will become" [7,9] are 
misleading. 

Before separately discussing the physical meaning of the different regions, it 
is interesting to ask ourselves why this hump and its associated correlated 
region has not been discussed in the past, giving rise to the aforementioned 
confusion. The only explanation we can offer is that the time series examined 
were simply too short, thus missing the timescales where correlations exist. It 
is instructive to notice that in Fig. 4 of Ref. [15], the power spectrum of the 
low drive case exhibits only regions A, B and C, while the spectrum for the 
high drive case barely extends to the beginning of region E. With that hmited 
amount of dynamical information it is reasonable that the existence of dy- 
namical correlations was mistaken for a consequence of the onset of avalanche 
overlapping. Proper credit, however, must indeed be given to the authors of 
Ref. [15] for their correct interpretation of the physical meaning of regions A, 
B and E. For that reason, we will just briefly discuss them in what follows, 
instead emphasizing the interpretation of regions C, D (the one that carries 
the signature of long time correlations) and F. 

4-1 Regions A and B (and A/B): Pulse Shapes 

Regions A and B at low drive and region A at high drive are due to the low 
level physics of the system that is manifested as the trapezoid pulse shape of 
individual avalanches (see Fig. 1). The values of the slopes and the location of 
the breakpoints of these regions can in fact be predicted by simply assuming 
a random superposition of trapezoids, in an analogous way to that which was 
originally done in Refs. [13,14] and noted for the high drive case in [15]. The 
reason why this approach works is that regions A and B do not depend on 
the character of the system dynamics. There are two high frequency regions 
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(A and B) in the sandpile spectrum as opposed to the single one derived in 
Refs. [13,14] because those studies assumed a rectangular pulse shape. Such 
pulses have a single spectral region while trapezoids have two. Simple tests 
show that replacing the trapezoids in the sandpile activity time series by 
either rectangular pulses or by spikes ('delta' functions) results in appropriate 
changes to the exponents in regions A and B but in no change to regions C, 
D and E. The R/S analysis exhibits only one region (A/B) that ends at a 
breakpoint that coincides with Tg ^. Therefore, regions A and B of the power 
spectrum and region A/B of the R/S analysis correspond to the same low level 
physics of the system. In fact, note that ~ 1 in region A/B, which implies 
that the value of the breakpoint Tb is related to the time series autocorrelation 
time and thus to the average avalanche duration. 

Therefore, it is clear that the value of the exponents found in regions A and B 
does not have a dynamical origin as has been pointed out previously [13,14,15]. 
This can be easily proved by shuffling the time series: keeping the trapezoids 
intact but placing them randomly in time. This reordering necessarily makes 
all dynamical correlations vanish. It is reflected in the spectra by the disap- 
pearance of regions D and E and the lack of change of regions A and B (see 
Fig. 7). This is an important statement of the lowest level physics of the sys- 
tem. It tells us that the high-frequency part of the spectra of systems with 
different low level physics and different fundamental pulse shapes may scale 
differently even if they exhibit the same low-frequency part of the spectrum 
because they share the same kind of dynamics. 



4-2 Region C: Quiet Times 

This is a region not previously identified in the running sandpile. In contrast 
to the vanishing drive limit, quiet times (periods with no activity) between 
avalanches appear in the running sandpile [5,6]. In essence, quiet times are a 
measure of the time that the system spends in an absorbing state. They are 
the source of the appearance of region C, which scales as /° in the power 
spectrum and as ~ 0.5 in the R/S analysis, both being signatures of an 
uncorrelated Poisson process. These exponents reflect the random character 
of the external drive that causes avalanches to be randomly triggered. Besides 
SOC models, quiet times have been studied in other physical systems [8,33]. 

Beyond the timescale of the largest single avalanche, Tb ~ L, there can be a 
period where the correlations are dominated by the random triggering. The 
width of region C is inversely related to the driving rate for a fixed system 
size. In fact, it disappears in the higher drive cases, after the average quiet 
time becomes on the order of the characteristic avalanche duration, Tb. Also, 
the lower the driving rate PqL, the fewer avalanches that occur in a given time 



13 




PX (effective driving rate) 



Fig. 8. H , (3, average quiet time and average duration versus five orders of magni- 
tude of effective driving rate. H stays constant with changing driving rate while /3 
increases as driving rate increases, saturating at around (3=1. 

period and the longer the system must 'wait' for enough avalanches to occur 
to correlate with each other. The cutoff, Tc, is a measure of this approximate 
minimum time and it scales inversely with the system time scale, ~ Pq^- 

Note that quiet times do not have a physical meaning in the zero drive limit 
of any sandpile, when the addition of sand is suspended during an active 
avalanche and is reinitiated immediately after it is completed. Their presence 
in the running sandpile does not change the dynamics, but allows the external 
drive to set a physical measure of time, Pq^, and to estimate the average quiet 
time [6]. The test shown in Fig. 7 justifies this statement: the spectrum and 
R/S analysis of the zero drive case are identical with those of the low drive 
case with quiet times removed. Note that a low drive case must be used for 
this exercise, so that there arc few or no overlapping avalanches. We emphasize 
this to introduce the next section, which shows that overlapping of events is 
not necessary to produce correlated dynamics. 

4-3 Region D: Dynamical Correlations 

Proving the existence of and understanding the physics of region D is the main 
point of this paper. For all cases in the low drive regime, regardless of system 
size, we find H ~ 0.8 in region D, indicating long time correlations. For the 
same timescale range, an region with < /9 < 1 is found in the power 
spectra, again a signature of correlations. While H remains constant, though, 
/3 changes with driving rate (see Fig. 8): it approaches 1 at the highest drives 
and decreases as drive strength is reduced. 

It is easy to show that the correlations must arise from interactions among 
distinct events and not from overlapping of events. The first clue is that this 
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Fig. 9. (Color online.) Power spectra and R/S for high drive with superpulses 
changed. Numbers shown are (3 for spectra and H for R/S. 

region exists on timescales far greater than the maximum duration of an in- 
dividual avalanche (Tb). Since the activity time series consists of a succession 
of avalanches separated by quiet times and there is no (or very little) over- 
lapping in the low drive regime, the correlations must be due to the specific 
order of events, chosen by the system dynamics in spite of the randomness 
of the drive. This is confirmed by the random shuffling of the activity series 
mentioned when discussing regions A and B (see Fig. 7). Regions C, D, E 
and F collapse into a single uncorrelated region for all time scales beyond Tb, 
where the shuffled series has a. spectrum down to the lowest frequencies 
and a Hurst exponent if ^ 0.5 to the longest time scales. 

A second confirmation that the correlations are established among the separate 
events is that quiet times have no role in the correlations. This is confirmed 
by the test already mentioned when discussing region C: removing all quiet 
times from the activity data while keeping the event order unchanged. As seen 
in Fig. 7, region C disappears: the H ~ 0.5 region is replaced by a if ~ 0.8 
region, which means that the H ~ 0.8 region D of the original data has moved 
to shorter timescales due to the removal of quiet times and the concomitant 
shortening of the series. But the correlations among events remain the same 
since neither their order nor their sizes have been changed. 

The same value of H ^ 0.8 is found for this region regardless of system 
size or driving rate, even in the vanishing drive limit in which SOC is strictly 
defined [17]. This simple fact suggests that the correlations present at the SOC 
critical point do persist almost without distortion at finite drives. In fact, this 
is still the case even at the highest drives, when avalanches do overlap and form 
superpulses. A superpulse is defined as the structure between successive quiet 
times in a high drive activity time series (Fig. 1(b)). Note that quiet times 
in the high drive cases are very small, on average. By randomly shuffling the 
superpulses of a high drive time series, we see that the correlations that lead 
to if ~ 0.8 are on time scales shorter than the average superpulse. Beyond 
these time scales, /? ~ and H ~ 0.5, signatures of uncorrelated data (Fig. 
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Fig. 10. Power spectra and (inset) R/ S analysis for fractional Gaussian noise of 
H = 0.8 with and without quiet times added. Numbers shown are [3 for spectra and 



Understanding the power spectrum of region D is not as simple as for the 
R/ S analysis. The slope of this region, /?, changes with driving rate and with 
system size. 1// regions are only obtained for high drives (or by removing the 
quiet times of low drive cases, which makes the system approach the zero drive 
limit as previously discussed) . This implies that a 1 // region is not a necessary 
signature of SOC dynamics at finite drives, since correlations are still present, 
as indicated by the Hurst exponent remaining constant at if ~ 0.8. In this 
sense, R/ S analysis is a more robust measure of the degree of correlations 
among events than is the power spectrum. 

The question, however, remains open about why (3 changes as region D is 
pushed towards different timescales, and why /? saturates at 1 and not another 
value. We do not know the answer, but we can say a couple of things. First, 
it appears that 1// is due more to a lack of quiet times and distinct pulse 
shapes than to overlapping, at least in the sandpile model. We showed in Fig. 
7 that removing all quiet times from between separate events of a low drive 
activity time series produces a spectrum with a region D where /3 = 1. In 
that case, there is no overlapping of events by design. We can also discard the 
need for hidden additional correlations being established dynamically, since 
the same change in (3 is observed when inserting randomly-distributed quiet 
times into an artificial time series of fractional Gaussian noise (fGn) [29] with 
an arbitrary Hurst exponent H = 0.8 (see Fig. 10). We inserted quiet times 
with an average approximately the same as the average of the original fGn 
series (~ 20). The spectrum and R/S analysis are shown for the original and 
modified series. The original fGn series follows /3 = 2H — 1, with H ~ 0.8 and 
13 ~ 0.6. When quiet times are added between each point of the original series 
a region of if ~ 0.5 appears up to a certain time lag. Beyond this time lag. 



H for R/S. 



9). 
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H ~ 0.8 as expected since the correlations among the data have not changed. 
In the spectrum, /5 ~ down to the inverse of the same time lag and then 
/3 ~ 0.43. This constant H and changing (3 is the same effect seen in the 
sandpile data. 

4-4 Region E: Global Discharge Region 

Region E is the discharge event region, aheady well-studied by Ref. [15] for 
high drive, and a consequence of the finite system size. The finite sandpile is al- 
ways being driven towards a globally critical slope. Once all or almost all cells 
have slope > Zcrit — ^f, a keystone toppling, usually near the bottom, will pro- 
duce a system-size avalanche or a rapid succession of smaller avalanches that 
removes enough sand so that the slopes at all sites are reduced to much less 
than critical (< Zcnt — Ni). Once such a large event occurs, it is unlikely that 
one will happen again for a long time. Hence, large events are anticorrclated 
and have the signatures of anticorrclated dynamics, /3 < in the spectrum 
and H < 0.5 in R/ S analysis. The same process drives both the low and high 
drive cases. 

4-5 Region F: Beyond Dynamics 

Beyond breakpoint TpT^, at the lowest frequencies, the spectrum is fiat, which 
we conjecture to be a refiection of the random drive of the system on the 
longest time scales. Note that all dynamics is on time scales much smaller 
than those in this region. Correlations are eventually erased by the system- 
wide events and shuffled by the random drive at time scales longer than Tp. 
Region F extends to infinitely low frequencies; there are no further dynamical 
regions below it. 



5 Relation between system size and driving rate 

We return to Figure 6 to discuss the effect of changing system size. It is rea- 
sonable to increase the level of forcing — practically, the amount of avalanche 
overlap — in the system by increasing the driving rate Pq. This means that 
more grains per time step will fall on the sandpile and triggering time will be 
reduced. It turns out that keeping Pq fixed while increasing system size L has 
the same effect because avalanches last longer in the larger system. In both 
of these cases, PqL has increased, causing trigger time to decrease. But since 
system size L has remained constant, the ratio of average avalanche duration 
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to trigger time has increased further, producing more overlap and stronger 
effective forcing. 

But the combination PqL does not by itself give enough information about 
the character of the drive as represented by the amount of avalanche overlap. 
Imagine increasing L while decreasing Pq so that PqL and, therefore, average 
trigger time remain constant. Though trigger time remains constant, the ratio 
of average avalanche duration to trigger time has again increased because of 
the increase in system size. This process, then, also effectively increases the 
driving rate in the larger system by producing more overlap among the indi- 
vidual avalanches that now last longer. This is true even though the original 
and the larger system both have the same PqL and, therefore, the same flux 
into the whole system and out of the bottom of the sandpile. 

In terms of the effective driving rate discussed in Section 2, Je is larger for the 
larger system when PqL is the same in both. This feature can be exploited to 
produce a plot qualitatively similar to Fig. 2. In that figure, the driving rate 

PqL of the sandpile increases from top to bottom for the spectra shown while 
system size L remains constant. If instead, one started with the spectrum 
of a small system at the top of the figure and increased system size while 
keeping PqL constant, then the same figure would be obtained. In both cases, 
PqL"^ increases from top to bottom. An application of this is that perhaps 
one wants to study a physical system that is suspected of being described by 
SOC and/or the running sandpile. If experimenters cannot change the system 
size (like, for instance, that of the magnetosphere) then perhaps periods of 
different observed driving rates could be used as a proxy for different system 
sizes. 

Since is the correct quantity to use to compare systems in the same drive 
regime but where the individual parameters are different, this has implications, 
for example, in investigations of finite size and/or multifractal scaling. As an 
illustration. Fig. 11 shows justification for comparing systems with the same 
effective driving rate, in this case with PqL'^ — 2.0. The spectra of systems 
with different driving rates PqL and system sizes L but with the same PqL^ 
can be rescaled to lie on top of each other. The spectra of systems with the 
same driving rate PqL but with different system size L cannot be rescaled to 
lie on top of each other. 

Finally, we show in Fig. 12 the PDFs of avalanche durations from two different 
system sizes that have been rescaled by a factor of L^" with a = 0.34. The 
rescaling has been done for a medium driving rate of PqL^ = 2.0 so that there 
is some avalanche overlapping. This, coupled with the rescaling of the spectra, 
shows that the running sandpile exhibits critical dynamics similar to that seen 
in the sandpile with vanishing drive [25]. 



18 




log frequency 

Fig. 11. (Color online.) Rescaling of power spectra of two different systems with 
PqL^ = 2.0, one with L = 200 and one with L = 2000, using the rescaling law 
yologiQ S{f) = g {xologio{f /xi)), where S{f) is the power spectrum and g{f) is a 
scaling function. Different sets of parameters [3;o,a;i,yo] are needed for each spec- 
trum so this is not the same as, for instance, the multifractal rescaling of avalanche 
size PDFs of [25], where the same value of each parameter is used for different 
system sizes. 
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Fig. 12. (Color online.) Rescaling of PDF of avalanche durations for two system 
sizes, L = 200 (solid line) and L = 2000 (dashed line), and a shared effective driving 
rate PqL'^ = 2.0. The PDFs have been rescaled by L° where a = 0.34. 

6 Conclusions 

We have analyzed the one dimensional directed running sandpile model for five 
orders of magnitude of driving rate and for almost three orders of magnitude 
of system size and have shown that the same dynamical correlations that are 
present in the SOC zero-drive limit persist (and produce nontrivial signatures 
in the power spectra andWe show through rescaling of power spectra and 
avalanche size PDFs that the running sandpile exhibits SOC in the same way 
that the vanishing drive sandpile does. R/S analysis of activity time series) 
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at finite drives. Our results warn that a time series from any system, for 
wliicli SOC may provide a plausible description of the dynamics, may be too 
short to see these correlation signatures and could thus be mistaken for a 
simple random time series. We show through rescaling of power spectra and 
avalanche size PDFs that the running sandpile exhibits SOC in the same way 
that the vanishing drive sandpile does.W 

The reason why these correlations persist away from the SOC critical point at 
zero drive can be understood when discussing the physical mechanism through 
which the sandpile establishes dynamical correlations. It is well known that 
memory in the sandpile is retained in the heights (and. therefore, the local 
gradients) of each cell. This is, for a physical system, the instantaneous system 
state. When an avalanche ends, the cells at the endpoints of its active zone 
are closer to critical and are more likely to be the initiation point of another 
avalanche. This is the basis for memory in the system and explains why the 
drive strength docs not affect dynamical correlations, except for pushing them 
to longer or shorter time scales: regardless of how long the intervening quiet 
time is between two events, they are correlated in the same way. In reference to 
the introduction, this is how an event can 'know' how large it will become — it 
can grow only as large as the size of the local neighborhood of marginal cells. 
Also, this is why an event of any size cannot simply occur at any time or place 
in the system. Only the time of the triggering of an event is random due to 
the random external forcing. Before that trigger, though, the size and location 
of the event is predetermined by the state of the local neighborhood, like an 
earthquake fault ready to slip. Therefore, the same memory mechanism that 
is the basis for SOC at zero drive is still present in the gradients regardless of 
how slowly or quickly sand is added, as long as the system is not overdriven. 

As we have shown, region D is the only region tliat reflects the underlying 
zero-drive SOC dynamics at finite drive: on the time scales in this region, the 
spectral and R/S signatures reflect only long time correlations and nothing 
about pulse shape, quiet times, random superpositions, overlapping of pulses 
or system size. In fact, in the thermodynamic limit in which the system size 
becomes increasingly large, we must simultaneously allow Pq to keep Je = 
PoL"^ constant, which pushes the system to the SOC critical point [17,18,19]. 
In that case, region D would extend to lower and lower frequencies. 

Finally, another interesting result derived from this work is that the relation 
P = 2H — 1 derived in Ref. [24] for fGn series and often accepted as the 
correct relationship between /3 and H regardless of the system, does not hold 
in general. While this fact is known from general theory [30,24] , the sandpile 
results provide a concrete example where the relation does not hold. This 
implies that the random drive Pq of the system for which f3 = 2H — 1 does 
hold is filtered by the sandpile so that (3 and H no longer have such a simple 
relation. Note that the changing /3 and constant H of the correlated region 
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D imply that the mechanism that produces 1// at the high drive produces 
non-1// at lower drive, so that '1// is not always 1//' and smaller values of 
f3 are not necessarily any less 'special' than f3 = 1. 
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